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ABSTRACT 


This  thesis  deals  with  the  enhancement  of  video  images 
degraded  by  turbid  water  viewing  conditions.  An  algorithm  by  Peli 
and  Lim  has  been  used  with  some  success  for  enhancement,  but  it 
was  found  to  accentuate  noise,  The  thesis  examines  a  combination 
of  the  Peli  and  Lim  algorithm  with  three  approaches  to 
enhancement. 

First,  a  Short  Space  Spectral  Subtraction  algorithm  which 
performs  the  restoration  in  the  density  domain,  using  an  estimate 
for  the  power  spectrum  of  the  given  data  set.  The  degraded  image 
is  divided  into  many  subimages  and  each  subimage  is  restored 
separately  and  then  combined. 

Next,  an  algorithm  for  Image  Enhancement  and  Noise  Filtering 
by  Use  of  Local  Statistics,  which  uses  the  assumption  that  the 
sample  mean  and  variance  of  a  pixel  is  equal  to  the  local  mean 
and  variance  of  all  pixels  within  a  fixed  range  surrounding  it. 

Finally,  a  median  filter  for  noise  reduction  ,  where  a  given 
pixel  of  a  degraded  image  is  replaced  by  the  median  of  the  pixel 
values  in  a  window  surrounding  it. 

Combination  of  the  algorithms  are  applied  to  degraded  images, 
and  the  results  are  compared  and  discussed,  in  each  case.  It  was 
found  that  noise  smoothing  can  be  achieved  with  the  spectral 
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I.  INTRODUCTION 


In  this  thesis,  some  algorithms  for  the  enhancement  of  video 
images  degraded  by  turbid  water  viewing  conditions  are 
implemented.  The  images  used  to  illustrate  the  processes,  were 
recorded  during  the  recovery  procedure  in  a  torpedo  testing  range. 

The  recover  operation  is  carried  on  undersea  with  the  aid  of 
digging  equipment  and  underwater  video  cameras  equipped  with 
strong  artificial  lights.  The  equipment  is  controlled  and  the 
operation  is  monitored  on  video  monitors  at  the  surface  onboard 
the  recovery  vessel.  The  recovery  equipment  in  its  attempt  to  dig 
out  the  torpedo,  stirs  up  sediment  wich  visual!/  obscures  the 
object  of  interest  and  impedes  the  operation.  To  date  no  special 
techniques  have  been  used  operationally  to  proce^  or  enhance  the 
video  image  before  display. 

When  it  is  of  interest  to  “improve”  an  imagp,  there  are  two 
broad  types  of  image  manipulation  processes  that  cover  all 
operations  performed  to  get  such  improvement:  image  restoration 
and  image  enhancement. 

The  goal  of  image  enhancement  is  to  process  a  degraded  image 
so  that  the  result  is  more  suitable  than  the  original  image  for  a 
specific  application  or  aids  the  human  analyst  in  the  extraction 
and  interpretation  of  pictorial  information  [Ref.  l].  Image 
restoration,  on  the  other  hand,  is  intended  to  bring  a  degraded 


image  back  to  an  ideal  degradation-free  image  as  closely  as  possible 
[Ref  2] .  This  makes  the  image  enhancement  problem  a 
subjective  one;  the  improvement  of  the  image  appearance  to  the 
human  viewer  is  highly  dependent  on  the  viewer  himself  and  on 
the  application.  What  is  “good”  for  one  person  or  application  is  not 
necessarily  “good”  for  another  one.  Consequently,  the  suitability  of 
the  processed  image  for  a  specific  application  makes  image 
enhancement  techniques  very  much  problem-oriented. 

The  degradation  process  can  result  from  causes  such  as,  an 
imperfect  photographic  process,  imperfect  display  devices,  poor 
contrast  due  to  environmental  conditions,  different  forms  of  noise 
(channel,  quantization,  salt  and  pepper)  and  others  [Ref.  ?].  In 
the  same  manner,  techniques  for  enhancement  are  dependent  of 
the  original  degradation  process.  Processes  yielding  satisfactory 
results  for  one  type  of  degradation  are  not  necessarily  suitable  for 
another. 

The  interest  of  this  thesis  is  on  examining  different  techniques 
to  enhance  images  that  are  degraded  by  environmental  conditions 
that  create  a  great  lack  of  contrast  in  the  recorded  images.  This 
work  is  a  continuation  of  previous  research  begun  on  the  subject, 
which  included  the  implementation  of  an  adaptive  filtering  [Ref.  4] 
algorithm  for  contrast  enhancement.  This  algorithm  was  applied  to 
a  set  of  underwater  images  recorded  in  turbid  water  viewing 
conditions  The  images  are  low  contrast  degraded  by  noisy 
background.  As  will  be  explained  in  the  next  chapter,  the 


algorithm  yielded  very  good  contrast  enhancement  but,  seemed  to 
accentuate  the  noise.  The  purpose  ot  this  work  is  to  investigate 
some  noise  smoothing  techniques  to  be  applied  to  the 
contrast-enhanced  images,  an  alternative  contrast  manipulation 
scheme,  and  to  compare  the  resulis  obtained  to  those  found  with 
the  algorithm  already  implemented.  Besides  the  images  included  in 
the  experimental  results  of  this  thesis,  many  more  images 
resulting  from  other  combinations  and  variations  of  the  processes 
and  different  settings  of  the  parameters  were  obtained.  We  have 


included  here  only  those  thought  to  be  most  representative  of  the 


results  obtained. 


A.  DESCRIPTION  OF  THE  ADAPTIVE  FILTERING  ALGORITHM 

An  adaptive  filtering  algorithm  for  image  enhancement  by 
Peli  and  Lim  [Ref.  4]  has  been  used  for  contrast  enhancement  of 
images.  The  images  are  typifyed  by  a  large  dynamic  range,  and 
recorded  over  a  medium,  characterized  by  a  much  smaller 
dynamic  range,  which  causes  those  regions  of  the  images  with 
very  high  or  very  low  luminance  to  be  poorly  represented  This 
algorithm  modifies  the  local  luminance  mean  of  an  image  and 
controls  the  local  contrast  as  a  function  of  the  local  luminance 
mean  of  the  image. 


f  (nl,n2) 
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Figure  2.1.  Adaptive  Filtering  for  Image  Enhancement 


In  Figure  2.1  a  block  diagram  of  the  Adaptive  Filtering 
algorithm  is  shown  In  the  figure,  f(ni,n2,'  denotes  the 
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unprocessed  digital  image,  and  f p (nt ,  n2) , which  denotes  the  local 
luminance  mean  of  f(n1,n2),  is  obtained  by  low-pass  filtering 
f(n1,n2).  The  sequence  fH(n1,n2)  which  denotes  the  local  contrast 
is  obtained  by  subtracting  fL(n1,n2)  from  i(n1,n2).  The  local 
contrast  is  modified  by  multiplying  fn(ni,n2)  with  k(fiJ,  a  scalar 
which  is  a  function  of  fL(ni,n2).  The  modified  contrast  is  denoted 
by  f'H(ni,n2).  The  specific  functional  form  of  k(fiJ  depends  on  the 
particular  application  under  consideration;  a  value  of  k(fjJ  >  1 
represents  local  contrast  increase  while  k(fiJ  <  1  represents  local 
contrast  decrease.  The  local  luminance  mean  is  modified  by  a  point 
non-linearity  and  the  modified  local  luminance  mean  is  denoted  by 
f'L(n1,n2).  The  specific  non-linearity  chosen  depends  on  the 

particular  application  under  consideration,  and  in  most  application 
problems  the  non-linearity  is  chosen  so  that  the  overall  dynamic 
range  of  the  resulting  image  is  approximately  the  same  as  the 
dynamic  range  of  the  recording  medium.  The  modified  local 
contrast  md  local  luminance  mean  are  then  combined  to  obtain 
g(n1,n2),  the  processed  image. 

B.  DISCUSSION  OF  RESULTS 

In  this  section  the  application  of  the  algorithm  described  in  the 
previous  section  to  the  images  shown  in  Figure  2.2  is  presented. 
Figure  2.2  shows  two  low-contrast  noise-degraded  images  of  512  X 
512  pixels  with  each  pixel  represented  by  8  bits.  Note  that  the 
small  dynamic  range  of  the  luminance  of  both  images  make 
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difficult  the  appreciation  of  subtle  details  in  them.  In  this  case,  it 
is  desired  to  increase  the  local  contrast  in  the  images  to  bring  back 
the  lost  details  To  achieve  this,  k(fL)  is  chosen  as  shown  in  Figure 
2.3  (a)  and  (b)  for  the  images  of  Figure  2.2(a)  and  (b), 

respectively.  For  the  case  of  Figure  2.2(a),  the  function  is  chosen 
to  increase  the  local  contrast  for  the  low  luminance  regions, 
maintaining  it  relatively  constant  for  medium  ranges  of  local 
luminance,  and  again  enhance  it,  for  the  brighter  regions  of  the 
degraded  image.  The  results  are  shown  in  Figure  2.4(a),  where  a 
relative  contrast  increase  can  be  observed.  This  is  especially 
evident  in  the  “fish”  shape  close  to  the  lower  right  corner  of  the 
image.  Also  note  the  increase  in  background  noise  resulting  in  the 
textured  background  observed  in  the  contrast-enhanced  image. 


Figure  2.3.  K(fJ  Function  for  Processing  Images  in  Figure  2.2 
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in  the  case  of  the  image  of  Figure  2.2(b),  it  is  desired  to 
retain  the  original  local  contrast  in  the  darker  regions,  but 


enhance  it  in  the  regions  corresponding  to  medium  and  high  local 
luminance  means.  This  explains  the  particular  shape  of  the  k(fj 
function  chosen  to  process  this  image.  In  the  results,  shown  in 
Figure  2. 4(b), some  details  difficult  to  be  distinguished  in  the 
un-processed  image,  are  now  more  evident.  Among  them  are  the 


particular  shape  of  the  object  near  the  center  of  the  image  and 
the  form  of  the  ocean  floor  below  the  left-most  object.  Again,  it 
can  be  observed  that  there  is  an  increase  in  the  background  noise 
throughout  the  contrast-enhanced  image. 
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III.  DESCRIPTION  OF  THE  ALGORITHMS 


In  this  chapter,  the  theorethical  foundations  and  procedures 
for  implementation  of  the  different  algorithms  used  in  the 
noise-filtering  and  contrast  enhancement  problem  are  presented. 
We  start  with  a  description  of  noise-smoothing  techniques  to  be 
used  as  post-processing  after  the  adaptive  filtering  algorithm.  Then 
an  alternative  method  for  the  contrast  enhancement  is  presented 
Finally,  the  use  of  median  filtering  for  noise-smoothing  is 
described 

A.  SHORT  SPACE  SPECTRAL  SUBTRACTION 

This  algorithm,  developed  by  J.  S.  Lim  [Ref. 5]  performs  the 
restoration  in  the  frequency  domain,  using  a  converging  solution 
for  the  power  spectrum  of  a  given  data  set  that  does  not  depend 
on  the  original  estimate.  The  procedure  is  based  on  two  basic 
assumptions.  First,  each  part  of  an  image  generally  differs 
sufficiently  from  other  parts  so  that  the  image  cannot  be  modeled 
by  a  stationary  random  process.  Second,  the  power  spectrum  of 
the  restored  image  is  estimated  by  the  spectral  subtraction  of  the 
additive  noise  spectra  from  the  degraded  image. 

1 .  Spectral  Subtraction 

Frequency  domain  techniques  for  image  enhancement  use 
iterative  procedures  for  estimating  the  image  power  spectral 


density,  continuing  until  certain  fidelity  criteria  are  met.  The 
density  so  obtained,  is  dependent  on  both  the  initial  estimate  for 
starting  the  iterative  procedure  and  the  stopping  criterion.  Some 
methods  converge  more  rapidly  than  others,  and  sometimes  first 
and  second  order  statistics  of  the  degrading  function  must  be 
known . 

Consider  a  model  of  image  degradation  as  presented  in 
Figure  3.1.  In  this  model  f(n1,ii2)  represents  a  digital  image, 
b(n1,n2)  represents  a  linear  space-invariant  point  spread  function, 
and  d(n1,n2)  represents  an  additive  noise  component.  Thus, 
g(n1,n2)  represents  a  blurred  version  of  f(ni,n2)  and  y(n1,n2)  is 
the  degraded  image. 


Figure  3.1.  Model  of  Image  Degradation 


An  image  restoration  system  is  shown  in  Figure  3  2.  It 
can  be  seen  that  the  goal  of  a  restoration  system  is  to  process  a 
degraded  image  y(n1,n2)  through  a  noise  reduction  system  to 
estimate  g(ni,n2)  which  in  turn  is  processed  by  a  deblurring 
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system  to  get  f(n1)n2),  an  estimate  of  the  original  blur-free  noise- 
free  image. 


Figure  3.2.  Image  Restoration  System 


The  purpose  of  this  algorithm  is  to  provide  a  model  for 
the  Noise  Reduction  System,  as  presented  in  the  Figure.  A  general 
model  for  the  restoration  filter  is  given  by 


w  -  { 


P  (0)  ,6)  ) 

9'  1  2 ' 


pg(w.«  )  ♦  «  •  Pd(“.“  ) 


(3.1) 


where  Pg(c«h,G)2)  and  Pd (&>i , (O2)  represent  the  power  spectra  of 
g(n1,n2)  and  d(n1,n2)  respectively,  a  and  P  are  constants.  If  P  is 
unity,  H((0i,(i)2)  corresponds  to  the  parametric  Wiener  filter,  and  if 
a  also  is  unity,  it  reduces  to  the  standard  Wiener  filter  [Ref.  1,2], 
Other  values  for  a  and  P  define  different  filtering  techniques,  such 
as  power  spectrum  filtering  and  geometrical  mean  filtering.  One 
commonly  used  approach  in  the  implementation  of  H(goi,co2)  is  an 

iterative  procedure  which  begins  with  an  initial  estimate  of 
Pg(co1,oo2)  and  then  iteratively  estimates  g(ni,n2)  and  Pg(u)i,<jo2) 
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until  a  converging  solution  or  a  desirable  performance  is  achieved 
This  iterative  procedure  may  be  computationally  undesirable. 

The  Spatial  Subtraction  approach  analytically  obtains  a 
converging  solution,  by  estimating  Pg(u>i,cj2)  from  g(n!.n2)  as 
(l/k.)* |G (ooi , 002) I2  and  using  a  value  for  |3  of  V2.  The  complete 
solution  for  this  problem,  is  given  by 


*  G(W  -W)  =  9  Y(W,a)  (3.3) 


for  |  Y(«1,«2)|2  i  a.k.P^U.tt) 


and  zero  otherwise 


i 


Here  G  (oi! ,  CO2)  represents  the  discrete  space  Fourier  Transform  of 
g(ni,n2),  Y , a>2)  is  the  discrete  space  Fourier  Transform  of 
y(n1,n2),  Pdfai,4^)  represents  the  power  spectrum  of  d(ni,n2), 
the  additive  noise  component  of  the  degradation  model,  and  k  is  a 
scaling  factor  that  normalizes  the  power  and  energy  spectral 
densities.  The  phase  of  g(n1.n2)  is  estimated  by  the  phase  of 
y(n1,n2)  and  the  transform  magnitude  of  g(n1;n2)  is  estimated  by 
a  particular  form  of  spectral  subtraction. 
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As  mentioned  before,  the  basic  assumption  is  that  each 
part  of  an  image  f(n1,n2)  generally  differs  sufficiently  from  other 

parts  so  that  it  cannot  be  modelled  by  a  stationary  random  field. 
Thus  a  short  space  implementation  seems  appropiate.  The 
technique  consists  of  dividing  the  original  image  into  many 
subimages  and  restoring  each  subimage  separately.  The  procedure 
involves  the  application  of  a  short  space  window  function 
wij(n1,n2)  to  overlapping  portions  of  the  degraded  image  y(n1(n2) 

so  that: 

y(ni'n2)'wiq(ni'n2)  =  {g(ni'n2)+d(ni'n2)}  *  vi,^ni'n2)  (3-4) 
or,  equivalently, 

yM(ni,n2)  =  (3.5) 

The  noise  reduction  system  is  then  applied  to  yij(n1,n2) 
to  recover  gyO'h,^).  The  overall  full-size  restored  image  g(ni,n2) 
is  given  by 


The  window  function  is  desired  to  be  a  smooth  function  to 
avoid  some  possible  discontinuities  that  may  appear  at  the 
subimage  boundaries  in  the  processed  image  and  must  satisfy: 

2K  2K 

2  2  wi  j(ni'n2>  =  1  (3-7) 

1-0  j-0 

for  all  rtj  of  interest 


The  condition  in  Equation  3.6  guarantees  that  an  image 
can  be  reconstructed  from  its  subimages.  To  ensure  the  smoothness 
of  the  window,  a  2-D  separable  triangular  window  similar  to  that 
shown  in  Figure  3.3  will  be  used.  Each  window  overlaps  its 
neighboring  window  by  half  the  window  kngth  in  each  dimension. 


B.  CONTRAST  ENHANCEMENT  BY  USE  OF  LOCAL  STATISTICS 


This  algorithm  presents  a  computational  technique  for  contrast 
enhancement  on  a  two-dimensional  image  array  based  on  the  local 
mean  and  variance.  The  algorithm,  as  developed  by  J.  Lee  [Ref 
6]  uses  the  basic  assumption  that  the  sample  mean  and  variance 
of  a  pixel  is  equal  to  the  local  mean  and  variance  of  all  pixels 
within  a  fixed  range  surrounding  it. 

Let  Xjy  be  the  brightness  of  a  pixel  (i,  j)  in  a  two  dimensional 

NXN  image.  The  local  mean  and  variance  are  calculated  over  a 
(2n  +l)x(2m  +l)  window.  The  local  mean  is  defined  as 


i+i 

_ 1 _ V 

“id  (2n+l)  '2n»i)  ^ 


(3.8) 


The  algorithm  is  designed  such  that  a  pixel  Xy  will  maintain 
its  local  mean,  and  yet  permit  its  variance  to  be  modif.ed  by  a 
constant  factor  times  its  original  variance  The  procedure  is  defined 


where  k,  the  gain,  is  the  ratio  of  new  local  standard  deviation  to 
the  original  standard  deviation.  This  approach  has  the 
computational  advantage  that  the  local  variance  Vjj  is  not 
required  and  only  the  local  mean  mjj  needs  to  be  computed  If 
we  have  k  >  1,  the  image  will  be  sharpened  as  if  acted  upon  by  a 
high-pass  filter.  If  0  <  k  <  1,  the  image  will  be  smoothed  as  if 
passed  through  a  low-pass  filter.  In  the  extreme  case, one  has  k  = 
0  and  x'j  j  is  equal  to  its  local  mean  mjj  . 


Figure  3.4.  Contrast  Enhancement  by  Use  of  Local  Statistics 


Figure  3.4  shows  a  block  diagram  of  the  algorithm.  The  local 
mean  is  evaluated  over  a  (2n+l)  (2m+l)  window  and  subtracted 
from  the  original  brightness  of  the  pixel,  xi(j.  The  difference  is 


then  multiplied  by  the  standard  deviation  ratio  to  enhance  the 
local  contrast.  Finally,  this  result  is  added  to  the  local  mean 
yielding  the  reconstructed  pixel  value  x'jj. 


C.  NOISE  FILTERING  BY  USE  OF  LOCAL  STATISTICS 

This  algorithm,  also  developed  by  J.Lee  [Ref.  6]  is  actually  a 
variation  of  the  previous  algorithm  adapted  for  the  noise  filtering 
problem.  In  this  algorithm,  the  a  priori  mean  (or  variance)  of 
the  estimated  image  is  calculated  as  the  difference  between  the 
mean  (or  variance)  of  the  noise  corrupted  image  and  the  mean 
(or  variance)  of  the  noise  itself.  Let  Zij  be  the  degraded  pixel  x*  j. 

Then  the  degraded  pixel  can  be  modelled  as  the  sum  of  the 
noise-free  pixel  and  a  white  noise  sequence.  That  is, 


z.  .  =  x,  .+  w.  . 
M  i/j  M 


(3.11) 


where  Wjj  is  the  white  random  sequence  with  E[wjj  ]=0.  Define 
the  a  priori  mean  and  variance  of  Xjj  to  be 


(3  12) 


°i.r  Ei<x1.rx‘i.)>  i=  Et(zi.rz"i.j) 


(3  13) 


respectively.  The  estimated  pixel  value  x^ j  is  computed  by 
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Figure  3.5.  Noise  Filtering  by  Use  of  Local  Statistics 

Figure  3.5  shows  a  block  diagram  for  the  implementation  of 
the  algorithm.  Starting  with  the  degraded  pixel  Zy.the  local  mean 

is  evaluated  and  subtracted  to  obtain  the  local  contrast.  Then, 
both  Qjj  and  ky  can  be  evaluated,  and  the  value  of  k  thus 

obtained  is  multiplied  by  the  local  contrast.  This  in  turn  is  added 
to  the  local  mean  to  get  the  restored  pixel  xy. 

One  important  parameter  is  the  size  of  the  window  over 
which  the  local  mean  is  estimated.  If  the  window  is  too  small,  the 


noise  filtering  algorithm  is  not  effective.  If  the  window  is  too  large, 
subtle  details  of  the  image  will  be  lost  in  the  filtering  process 

D.  MEDIAN  FILTERING 

Median  Filtering  is  a  non-linear  signal  processing  technique 
originally  developed  by  J.Tukey  [Ref.  7],  that  is  useful  for  noise 
suppression  in  images.  In  one-dimensional  form  the  median  filter 
consists  of  a  sliding  window  encompassing  an  odd  number  of  pixels. 
The  center  pixel  in  the  window  is  replaced  by  the  median  of  the 
pixels  in  the  window.  The  median  of  a  discrete  sequence  a1>a2)  .  .  ., 

aN  for  M  odd  is  that  member  of  the  sequence  for  which  (M-l)/2 
elements  are  smaller  or  equal  in  value,  and  (M-l)/2  elements  are 
larger  or  equal  in  value  [Ref.  2],  For  example,  if  the  ordered 
values  of  the  pixels  within  a  window  are  80,  90,  200,  110,  120, 
the  center  pixel  would  be  replaced  by  the  value  110,  which  is  the 
middle  point  of  the  sorted  sequence  80,90,110,120,200. 

The  concept  of  median  filter  is  easily  extended  to  two 
dimensions  by  utilizing  a  two-dimensional  window  of  some  desired 
shape  such  as  a  rectangle  or  a  discrete  approximation  to  a  circle. 
A  two-dimensional  LxL  median  filter  will  provide  a  greater  degree 
of  noise  suppression  than  sequential  horizontal  and  vertical 
processing  with  LXl  median  filters,  but  two-dimensional  processing 
also  results  in  greater  signal  suppression. 

The  median  filter  has  been  found  to  be  more  effective  than  a 
linear  filter  for  smoothing  images  with  spiky  noise  degradations 


because  of  extrema  rejection  by  the  median  [Ref.  8] .  Furthermore, 
the  median  filter  preserves  monotonic  step  edges,  that  is,  it  does 
not  blur  sharp  edges  as  a  linear  low-pass  filter  would. 

Another  interesting  aspect  of  median  filtering  that  has  been 
studied  [Ref.  9],  is  that  of  the  convergence  of  a  filtered  image  to 
what  has  been  called  a  root  signal  .  A  root  is  a  signal  which  is 
invariant  under  filtering  by  a  particular  median  filter.  The 
hiuque  consists  in  making  succesive  passes  of  a  noise-degraded 
image  thru  a  median  filter  until  an  image  corresponding  to  a  root 
signer  is  achieved.  Definition  of  a  root  signal  follows,  where  a 
window  width  of  2N+1  is  used. This  definition  uses  the  following 
ideas 

1)  A  constant  neighborhood  is  a  region  of  at  least  N+l  consecutive 
identically  valued  samples. 

2)  An  edge  is  an  increasing  or  decreasing  sequence  of  samples 
which  is  immediately  prece  ’  ^  and  followed  by  constant 
neighborhoods.  An  edge  cannot  contain  any  constant 
neighborhood. 

3)  An  impulse  is  a  sequence  of  at  most  N  consecutive  samples 
whose  values  are  different  from  those  of  the  two  surrounding 
regions;  the  two  surrounding  regions  are  identically  valued 
constant  neighborhoods . 

A  signal  is  a  root  if  and  only  if  it  contains  only  edges 
alternating  with  constant  neighborhoods  [Ref.  10],  The  median 
filtering  preserves  edges  and  constant  neighborhoods  but  eliminates 
impulses.  The  technique  is  based  on  the  idea  that  passing  a  noise 
corrupted  root  signal  through  the  median  filter  a  sufficient  number 


of  times  will  produce  another  root  signal,  which  is  usually  “close” 
to  the  original  root  signal.  In  reference  [Ref.  10],  the  concept  of 
“closeness”  and  a  measure  of  the  number  of  runs  to  be  used  are 
calculated  as  a  function  of  the  image  and  window  sizes. 

In  this  thesis,  a  median  filter  included  in  the  Spider  library 
[Ref.  11]  is  used,  which  is  the  implementation  of  the  algorithm 
suggested  by  Huang,  Yang  and  Tang  [Ref.  12].  It  consists 
primarily  of  an  efficient  way  of  updating  the  histogram  each  time 
that  a  given  pixel  is  replaced  by  the  median  within  a  window.  The 
algorithm  is  described  briefly  for  a  3X3  window: 

Step  1:  Obtain  the  histogram  in  the  first  window  [Figure  3.6 

(a)]  and  find  the  median  MDN.  Next, count  the  pixels 
with  values  smaller  than  the  median  of  the  window. 
LTMDN  is  that  number. 

Step  2:  Shift  the  window  right  by  one  pixel  [Figure  3.6(b)] 

and  update  the  histogram  and  LTMDN.  First  decre¬ 
ment  the  counted  values  of  the  histogram  equivalent 
to  gray  level  values  {gJ (a) gl (d) ,  gl (g)  )  of  pixels  a,  d 
and  g.  That  is: 

HIST [gl(a)  ]  =  HIST [gl(a)  ]-l 
HIST \gl(d)  ]  =  }UST[gl(dJ  ]-l 

hist[^;  ]  =  hist[^;  ]-i 

LTMDN  is  updated  as  follows: 

if  gl(a)  <  MDN  ,  LTMDN  =  LTMDN- 1 
IF  gl(d)  <  MDN  ,  LTMDN  =  LTMDN- 1 


if  g’(g)  <  MDN  ,  LTMDN  =  LTMDN-1 
In  the  same  way  for  pixels  c',f',  and  i‘ , 
HIST [gl(c)  ]  =  HIST [gl(c)  ]  +  l 
HIST [gJ(f)  ]  =  HIST [gl(f)  ]+l 
HIST [gl(l)  ]  =  HIST [gj(l')  ]+l 
IF  gl(c')  <  MDN  ,  LTMDN  =  LTMDN +1 

IF  gl(f')  <  MDN  ,  LTMDN  =  LTMDN+1 

IF  gl(i')  <  MDN  ,  LTMDN  LTMDN-H 


Figure  3.6.  Window  for  Median  Filtering 


Step  3:  Update  MDN  in  the  previous  window,  and  obtain  a 

median  in  the  new  window.  Let  ITH=  (number  of 
pixels  in  the  window)/2 
-  If  LTMDN  >  ITH,  then 

LTMDN  =  LTMDN  -  hist[MDN] 


MDN  =  MDN-1 

Repeat  this  until  LTMDN  <  1TH  is  obtained. 

-  If  LTMDN  <;  ITH,  then 
MDN  =  MDN+1 

LTMDN  =  LTMDN  +  hist[MDN] 

Repeat  this  until  LTMDN  +  hist[MDN+1]  >  ITH  is 
obtained.  This  MDN  is  the  median  for  the  current 
window  and  substitute  it  in  the  output  image. 

End  when  one  line  is  finished.  Go  to  Step  2. 


In  this  chapter,  experimental  results  for  different  combinations 
of  the  algorithms  described  are  presented  and  discussed.  Basically, 
the  adaptive  filtering  algorithm  [Ref.  4]  was  used  for  contrast 
enhancement,  as  described  in  Chapter  II.  As  previously  stated, 
this  algorithm  yielded  very  good  contrast  enhancement,  but  also 
tended  to  accentuate  the  noise.  Thus,  methods  for  noise  reduction, 
such  as  short  space  spatial  subtraction  [Ref.  5],  noise  filtering  by 
use  of  local  statistics  [Ref.  6],  and  median  filtering  are  used,  as 
post-processors  to  the  contrast  enhancement  operation.  The  results 
are  shown  in  next  sections.  Another  technique  for  contrast 
manipulation, contrast  enhancement  by  use  of  local  statistics  [Ref. 
6]  was  also  used,  as  an  alternative  to  the  adaptive  filtering 
algorithm.  These  results  are  also  given  and  compared  to  those 
previously  obtained. 

A.  PROCESS  TYPE  1 

This  section  shows  a  combination  of  the  adaptive  filtering  and 
short  space  spatial  subtraction  algorithms.  In  Figure  2.2  two 
degraded  images  are  shown.  The  images  are  characterized  by  poor 
contrast  due  to  the  turbid  water  viewing  conditions.  These  images 
were  processed  for  contrast  enhancement  using  the  adaptive 


filtering  technique,  but  the  noisy  background  was  heavily 
accentuated  as  shown  in  Figure  2.4.  Then  the  spectral  subtraction 
algorithm  was  used  to  smooth  the  noise  resulting  in  the  images 
shown  in  Figure  4.1.  Inspection  of  these  images  shows  that  this 
particular  algorithm  was  succesful  in  smoothing  the  background 
noise  present  in  the  contrast-enhanced  images,  and,  at  the  same 
time,  most  of  the  information  is  preserved.  Particularly  note,  in 
Figure  4.1(a),  how  details  in  the  beam-shaped  object  in  the  center 
of  the  image  are  brought  back,  as  well  as  those  of  the  “fish”  in 
the  lower  right  corner  of  the  image.  The  same  observation  can  be 
made  about  Figure  4.1(b),  in  which  the  noise  reduction  was 
achieved, at  the  cost  of  some  signal  degradation,  in  this  case. 

This  algorithm  requires  the  power  spectral  density  of  the  noise 
to  be  estimated.  In  this  case,  this  was  accomplished  by  taking  a 
portion  of  the  noisy  background,  without  any  objects  and  using  it 
to  estimate  the  spectrum.  Figure  4.2  shows  a  representation  of  the 
estimated  spectrum  of  the  noise  density  thus  obtained.  In  the 
figure,  the  magnitude  of  the  power  spectrum  is  represented  as  the 
intensity  of  a  2-D  image,  each  of  the  image  dimensions  being  the 
corresponding  spatial  frequency  coordinates. 

For  the  short  space  implementation,  it  is  required  to  process 
the  degraded  image  by  dividing  it  into  subimages  and  processing 
each  subimage  separately.  In  dividing  the  degraded  image,  the  size 
of  the  subimage  must  be  such  that  the  image  luminance  in  it  can 
be  approximated  by  a  stationary  random  process.  In  addition,  the 


Figure  4.1.  Images  in  Figure  2.4  Processed  with  the  Short 
Space  Spatial  Subtraction  Algorithm 


Figure  4.2.  Power  Spectrum  of  the  Noise  in  Figure  2.4  (a) 

window  function  must  be  a  smooth  function  in  order  to  avoid 
possible  discontinuities  that  may  appear  at  the  subimage  boundaries 
in  the  processed  image.  Thus,  two  types  of  smooth  functions  were 
used-  separable  2-D  triangular  and  Hanning  windows-of  size  32x52 
pixels,  overlapped  with  its  neighboring  window  by  half  the  window 
duration  in  each  dimension.  It  was  noted  that  there  was  no 
appreciable  difference  in  using  either  window  -triangular  or 
Hanning. 

The  parameter  k,  appearing  in  Equation  3.1  is  a  scaling  factor 
that  normalizes  the  power  spectral  densities,  and  was  evaluated  as 

OO  00 

k  -  2  2  MW  (11) 
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where  wO*,^)  is  the  window  function  used. 


In  this  section,  a  combination  of  the  adaphve  filtering  for 
Image  enhancement  and  noise  filtering  by  use  of  local  statistics  is 
presented.  The  images  in  Figure  2.4  were  processed  with  the  noise 
filtering  algorithm,  using  different  values  for  the  noise  variance 
Based  in  the  data  used  for  generating  the  noise  power  spectrum 
shown  in  Figure  4.1,  a  noise  variance  of  60  was  estimated,  for  the 

image  in  Figure  2.4(a),  For  the  image  in  Figure  2.4(b),  a  noise 

variance  of  45  was  estimated.  The  algorithm  was  found  to  be  quite 
sensitive  to  the  size  of  the  window  used  to  evaluate  the  local 

statistics,  the  a  posteriori  mean  and  variance.  Window  sizes  of 
8X8,  16X16,  and  24X24  pixels  were  used  without  noticing  any 
appreciable  improvement  in  the  processed  image.  Figure  4.3  shows 
the  results  of  processing  the  images  in  Figure  2.4  using  a  window 
size  of  32X32  and  the  noise  variances  mentioned  above.  Note  that 
a  relative  reduction  in  the  background  noise  was  achieved,  with  the 
processed  image  resembling  the  effect  of  low-pass  filtering  in  the 
background,  but  without  blurring  the  edges  of  the  objects  in  the 
images.  When  larger  sizes  for  the  window  were  used,  some 

distortion  in  the  images  begun  to  appear, 

C.  PROCESS-TYPE  5 

In  this  section,  the  result  of  applying  a  median  filter  to  a 
contrast-enhanced  image,  using  the  adaptive  filtering  algorithm,  is 
presented  The  image  in  Figure  2.4(a)  was  processed  with  the 


median  filter  available  in  the  Spider  Library,  using  a  window  size 
of  5X5,  and  making  succesive  iterations  of  the  processed  image, 
until  appreciable  image  degradation  start  to  show  up.  The  image 
obtained  after  8  iterations  is  shown  in  Figure  4.4.  It  can  be 
observed  that  some  noise  smoothing  was  obtained.  Even  when  the 
resulting  image  is  not  noise-free,  the  remaining  noise  is  of  much 
lower  spatial  frequency  than  that  of  the  image  before  the 
median-filtering,  and  thus  seems  less  objectionable. 

Detailed  characteristics  of  the  signal  are  preserved  and  some 
details  have  actually  been  enhanced.  This  is  true  in  the  case  of  the 
fish  in  the  lower  right  corner  of  the  image.  A  comparison  with  the 
use  of  3X3  window  for  the  filter  operation  led  to  the  conclusion 
that  the  smaller  window  size  gives  slightly  more  fidelity  in  the 
signal  preservation.  However,  since  convergence  for  the  3X3 
window  is  much  slower,  this  choice  is  undesirable.  A  reduction  in 
the  dynamic  range  of  the  median-filtered  image  can  also  be 
observed,  which  presents  an  effect  similar  to  that  of  low-pass 
filtering  in  the  background  region  of  the  image.  However,  this 
procedure  does  not  have  the  undesirable  characteristic  tendency  of 
the  low-pass  process  to  blur  the  edges  of  the  signal. 


Figure  4.4.  Median  Filtering  of  Image  in  Figure  2.4(a) 


D.  PROCESS  TYPE  4 

In  this  section  the  application  of  the  contrast  enhancement  by 
use  of  local  statistics  algorithm  is  presented.  In  Figures  4.5  and 
4.6  four  low-contrast  noise-degraded  images  are  shown.  All  of 
them  are  512  X  512  pixels  in  size  and  each  pixel  is  represented  by 
8  bits.  The  images  in  Figures  4.5  (a),(b)  and  4.6(a)  are 

characterized  by  a  relatively  dark  background  with  a  diffused 
effect,  around  the  objects  in  the  pictures  resulting  in  low  contrast. 
The  image  in  Figure  4.6(b)  has  an  almost  uniform  brightness 
;  throughout  the  picture,  which  makes  it  difficult  to  distinguish 

i 

between  the  fish  and  beam  shapes  from  the  uniform  background. 

Figure  4.7(a)  corresponds  to  the  image  in  Figure  4.5(a) 
processed  with  the  contrast  enhancement  algorithm,  using  a  value 

I 
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of  4  for  the  gain  k  of  Equation  3.10.  (Recall  that  k  is  the  ratio  of 
the  iocal  standard  deviation  of  the  processed  image  to  the  original 
standard  deviation) .  Note  the  enhancement  in  contrast  achieved, 
which  is  particularly  evident  in  the  cable  connecting  the  two 
objects  in  the  lower  part  of  the  image. 

The  image  in  Figure  4.5(b)  was  processed  using  a  va'ue  of 
k=6,  with  results  shown  in  Figure  4.7(b).  The  higher  value  was 
used  to  try  to  enhance  the  fish  shape  in  the  lower  right  part  of 
the  image.  This  enhancement  was  achieved  to  some  extent. 

Figure  4.8(a)  shows  the  result  of  processing  „he  image  in 
Figure  4.6(a).  In  this  case,  a  value  for  k  of  4  was  again  used, 
producing  an  image  that  is  a  much  crisper  version  of  the 
corresponding  degraded  image.  Also  note  that  details,  such  as  the 
particular  shape  of  the  object  in  the  picture,  are  enhanced.  The 
noise  in  the  background  is  not  greater  than  that  in  the  original, 
which  means  that  the  desired  result  of  contrast  enhancement  was 
achieved  without  a  significant  noise  increase. 

The  results  of  processing  the  image  in  Figure  4.6(b)  are  shown 
in  Figure  4.8(b).  In  this  case  a  value  of  k=5  was  used,  and  the 
contrast  of  the  two  objects  in  the  picture  are  increased  with 
practically  no  increase  in  the  noise. 
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Figure  4.6.  Original  Degraded  Images 
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Restoration  of  images  degraded  by  turbid  water  was  considered 
in  this  thesis.  Previous  work  on  the  subject  included  the  use  of  an 
adaptive  filtering  algorithm  for  enhancing  the  local  contrast.  This 
algorithm  was  applied  to  images  characterized  by  poor  contrast  due 
to  the  turbid  water  viewing  conditions.  The  technique  yielded  good 
contrast  enhancement,  but  tended  to  accentuate  the  noise. 

In  our  work  some  techniques  for  noise  smoothing  were  used  as 
post-processors  to  the  contrast  enhancement  operation. 

The  specific  techniques  used  were  short  space  spectral 
subtraction,  noise  filtering  by  use  of  local  statistics,  and  median 
filtering.  The  spectral  subtraction  algorithm  yields  very  good  noise 
reduction,  but  tends  to  introduce  some  signal  degradation.  The 
filtering  technique  that  uses  local  statistics  -local  mean  and 
variance-  produced  noise  smoothing  in  the  background.  This 
technique  was  found  to  be  very  sensitive  to  the  size  of  the  window 
in  which  the  local  statistics  were  evaluated.  Median  Filtering  did 
not  eliminate  the  noise  completely,  but  the  remaining  noise  was  of 
much  lower  spatial  frequency  than  the  original  noise  and  thus 
seemed  less  objectionable.  Better  noise  filtering  was  achieved  when 
the  degraded  image  was  median-filtered  iteratively  until  some 
signal  degradation  began  to  appear. 


As  an  alternative  to  the  adaptive  filtering  algorithm,  an 
algorithm  for  contrast  enhancement  by  use  of  local  statistics  was 
used.  This  particular  technique  yielded  very  good  results  without 
significant  noise  increase. 

All  processes  implemented  in  this  thesis  proved  to  be  quite 
computationally  intensive.  This  was  especially  true  for  the 
techniques  that  require  the  evaluation  of  local  statistics.  The 
ultimate  goal  for  the  problem  addressed  in  this  and  previous  work, 
is  to  perform  the  enhancement  in  real  time.  This  suggests  a 
natural  continuation  for  the  work  done  so  far  on  this  topic,  in 
order  to  implement  the  algorithms  for  enhancement  in  real  time. 
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APPENDIX 


COMPUTER  PROGRAMS 


PROGRAM  SPATSUB 

C 

c 

C  THIS  PROGRAM  IS  THE  IMPLEMENTATION  Of  THE  SHORT  SPACE 

C  SPECTRAL  SUBTRACTION  ALGORITHM 

C 
C 
C 

BYTE  BYTEIMG(512,5I2) 

CHARACTER"  1  COEN 

CHARACTER- 40  INFILENAME, OUTFILENAME 

INTEGER-4  IMGBUFF(5 1 2.5 1 2).  INTARRAY(5I2,5I2) 


C  *  -  »  INPUT  FILENAMES  :  INPUT  &  OUTPUT  -  *  - 

TYPE  10 

10  FORMAT  C  INPUT  FILENAME  =>  *.$) 

ACCEPT  15.  INFILENAME 
15  FORMAT  (A40) 

C  ACCEPT  67.  IWY 

67  FORMAT  (12) 

TYPE  68 

68  FORMAT  f  COLUMN  SIZE  OF  WINDOW  «>  '.$) 
IWX  =16 

C  ACCEPT  69.  IWX 

69  FORMAT  (12) 


I  ERR  =  0 

-«  -  INPUT  ORIGINAL  IMAGE  *"« 

CALL  INPUTIMG  (INFILENAME, BYTEIMG,  IMGX,  IMGY) 

-  -  -  CHANGE  IMAGE  FROM  BYTE  TYPE  TO  INTEGER  TYPE  -  -  » 

CALL  BYTE-TOJNTEGER  (BYTEIMG,  IMGBUFF,  IMGX.  IMGY,  IERR) 
IF  (IERRR  EQ.  1)  GO  TO  100 

•  -  -  PERFORM  THE  SPATIAL  SUBSTRACTION  *  -  * 


CALL  CORTO  (IMGBUFF,  INTARRAV,  IMGX,  IMGY,  IWY,  IERR) 


IF  (IERR  .EQ.  I  )  GO  TO  100 

--*  CHANGE  IMAGE  FROM  INTEGER  TYPE  TO  BYTE  TYPE  --* 

CALL  INTEGER_TO_BYTE  (INTARRAY,  BYTEIMG,  IMGX,  IMGY, IERR) 

IF  (IERR  .EQ.  1)  GO  TO  100 

*  *  *  OUTPUT  THE  PROCESSED  IMAGE  *  *  * 

CALL  OUTPUTING  (OUTFILENAME,  BYIEIMG,  IMGX,  IMGY) 

STOP 

CONTINUE 

TYPE  !!  FORCED  TO  EXIT  !!  1 
STOP 

END 

SUBROUTINE  CORTO  (ENTRAOA,  SALIDA.iri6X.IM6Y.IWY.IERR) 

THIS  SUBROUTINE  PERFORMS  THE  SPATIAL  SUBTRACTION  ALGORITHN 
USING  THE  SHORT  SPACE  IMPLEMENTATION  TECHNIQUE  FOR  IMAGE 
RESTORATION 


INTEGER-  A  ENTRADA  (IMGX.IMGY),  SALIDA(IMGX.IMGY),  INTEC32.32) 

REAL" 4  AR(32,32),AI(32,32),BR(32,32),BI(32,32), 

STM),  CT(64),  LBRM),  BRR(32.32),  WIND0W(32,32). 
RENTRADAC5 1 2,5 1 2),  IMGWIN(32,32),  IARRAY(32,32),  MAG(32,32), 
PHAS(32,32),PSD(32,32) 

DO  I  -  1  ,IMGX 
DO  J  =  1  .IMGY 

RENTRADA  (l,J)  -  FLOATJ(ENTRADA(l,J)) 

END  DO 
END  DO 

TYPE  -.INTEGER  IMAGE, CONVERTED  TO  REAL  DATA' 

CALL  MCHECK(RENTRADA.IMGX) 

-»»  BUILDUP  THE  WINDOW  FUNCTION  «** 

CALL  BLTWINDOW  (16.WIND0W) 

TYPE  -  .WINDOW  FUNCTION  BU'LT' 


CALL  MCHECK(WIND0W,32) 


r\ 


M 


AK  =0.0 


00  I  -  1 ,32 
00  J  -1,32 


AK  =  AK+(WIND0W(!,J)**2 


END  DO 
END  DO 


TYPE  * ,AK">...  \AK 


*  *  *  INITIALIZE  OUTPUT  ARRAY  *  *  * 


DOL  =  1  .IMG 
DO  M  -  I  ,IM6Y 

SALIDA  (L.M)  -  0 

END  DO 
END  DO 


TYPE  * , 'OUTPUT  ARRAY  INITIALIZED' 


CALL  MCHECK(SALIDA,IM6X) 


*  *  *  ESTIMATE  POWER  SPECTRUM  DENSITY  OF  THE  NOISE  *  *  * 


CALL  ESTINO!SE(RENTRADA,PSD) 


TYPE  * /POWER  SPECTRUM  DENSITY  Of  THE  NOISE  ESTIMATED' 


CALL  MCHECK(PSD,32) 


MAIN  PROCESS  SECTION  «< 


INSERT  LIMITS  HERE 


DO  I  -  1.31 
DO  J  -  1,31 

LL  -  (0-0*  16)+ 1 
MM  -  ((J-1)*  16)+ 1 
DOL  -  1,32 
DOM  -  1,32 

IA  - LL+L-1 
B  -  MM+M-l 

IFOA  ,6E.  16  .AND  IA  IE.  112  .AND. 

IB  ,6E.  16  .AND.  IB  IE.  112)  THEN 

IM6WIN(L,M)*RENTRADA(IA,IB) 

IARRAY(L,M)-WIND0W(L,M)*IM6WIN(L,M) 


ELSE  IFOA  .GE.16  .AND.  I A  ,LE.  112)  THEN 

IARRAY(L,M>*  0.5*RENTRADA0A,IB) 
ELSE  IFOB  ,6E.  16  .AND.  IB  IE.  112)  THEN 

IARRAY(L,MW).5*RENTRADA0A,IB) 


ELSE 

IARRAY  (L,M)=RENTRADA(IA,IB) 

END  IF 

END  DO 
END  DO 

TYPE  *  .DATA  ARRAY  HAS  BEEN  WINDOWED' 

CALL  MCHECK  (IARRAY.32) 

*  *  *  TRANSFORM  DATA  ARRAY  INTO  COMPLEX  FORMAT  *  "  * 

DO  L=1,32 
DO  M=  1 .32 

AKL.M)  =  0. 

AR(L.M)  =  lARRAY(L.M) 

END  DO 
END  DO 


TYPE", 'CHECK  REAL  PART' 

CALL  MCHECK  (AR.32) 

TYPE  ".'CHECK  IMAGINARY  PART  (ZEROS)' 

CALL  MCHECK  (AI.32) 

"  "  *  EVALUATE  2-D  FFT  OF  DATA  IMAGE  "  "  " 

CALL  FFTS2(AR,AI.BR,B),32>32,ST,CT,LBR.64.2.JERR) 

CALL  MCHECK  (8R.32) 

DO  IX  =  1.32 
DO  IY  =  1.32 

BRRdX.lYKABS  (CMPLX(BR(IX.IY).BI(IX.IY))) 

END  DO 
END  DO 

TYPE  ".TFT  DONE, CHECK  MAGNITUDE' 

CALL  MCHECK  (BRR.32) 

CALL  DISPECT  (BRR) 

"  "  "PERFORFM  SPATIAL  SUBSTRACT10N  "  "  " 

ALFA  =  1 .5 
PSD  =  40.0 
DO  IX  =  1.32 
DO  IY  =  1,32 

BR(IX.IY)  =  ((BRR(IX.IY))*  *2)  -  ((ALFA/AK)  »  ((PSD(IX.IY)  " 


IF  (BR(IX.IY)  LT.  .0.)  THEN 
MAG  (IX, IY)  =  0 

ELSE 

MAG(IX.IY)  =  SQRT(BR(IX.IY)) 

END  IF 

END  DO 
END  DO 

TYPE  *,' SPATIAL  SUB T. DONE  CHECK  RESULT' 

CALL  fiCHECK  (MAG) 

"  *  *  EVAL ,  MAGNITUDE  OF  INV,  FFT  *  *  * 

Dr  IX=  1,32 
DO  IY-  1,32 

BR(IX.IY)  =  MAGOX.I  Y )"  COS(PHAS(IX,l  Y ) 

Bl  (IX, IY)  =  MAG(IX.IY)* SIN(PHAS(IX,IY) 

END  DO 
END  DO 

CALL  FFTS2(BR,BI,AR.Ai,32,32,ST.CT,LBR.64.-2.JERR) 

TYPE*. INV.  FFT  DONE,  CHECK  REAL  PART' 

CALL  MCHECK(AR) 

TYPE*  .CHECK  IMAGINARY  PART  (ZEROS)  ' 

CALL  MCHECK(AI) 

*  *  *  BACK  TO  INTEGER  FORMAT  *  “  * 

DO  IX  =  1,32 
DO  IY  =  1,32 

INTE(IX.IY)  =  JINT(AR(IX.IY)) 

END  DO 
END  DO 

TYPE*, 'INTEGER  CONVERSION  DONE, CHECK  IT' 

CALL  MCHECK  (INTE) 

*  *  *  BUILD  UP  THE  OUTPUT  ARRAY  *  *  * 

DO  IX  =  1,32 
DO  IY  =  1 .32 

LL  =((I-U*  16}+ 1 
MM((J-1)*  16H-1 
IXX  =LL+IX-1 
IYY  =MM+IY-1 

SALIDA  (1XX.IYY)  ■  SALIDA(IXX,IYY)+INTE(IX,IY) 
TYPE*, 'IXX  AND  IYY  =>',IXX,IYY 
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m 


DO 

DO 


END  DO  !!  END  DO  J 
END  DO  !!  END  DO  I 

100  TYPE  M1AIN  PROCESS  DONE' 

C  CALL  TRUNCATE(SALIDA,IMGS,IMGY) 

CALL  SCALE_ADJUST(SALIDA,IMGX,IMGY, 0.0, 255.0) 
TYPE  *  ,’SPATIAL  PROCESS  DONE' 

RETURN 

END 


SUBROUTINE  BLTWIND0W(IWY.NIN2) 

REALM  NK32),  NIN2C32.32),  N2(32,32) 


L  =0 

PI  .  4-ATAN  (1.) 

DO  I  +  1 ,32 

II  =1-1 

N 1  (I )  *  0,5  -  (0.5** C0S(2W PI "  11/3 1 )) 
N2(l)  =  NKI) 

TYPE  *,  NKI),  N2(l) 

IF  (I  .ST.  16)  THEN 
L  =L+1 
NKI)  =  I  -  (2*L) 

ELSE 


NKI)  = 
END  IF 


END  DO 


L  =0 


C 

C 

C 

C 

C 

C 

C 


DO  1=1,  32 
DO  J  =  1 , 32 

IF  (J  .GT.  16)  THEN 
L  =L  +1 

NIN2(U)  -  NKI)  * 
ELSE 

NIN20.J)  =  NKI)  * 

END  IF 

NIN20.J)  =  N1(I)*N2(J) 
TYPE  \  NIN20.J) 

END  DO 
L  =0 


(J-2-1) 

J 
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RETURN 

END 


SUBROUTINE  ESTINOISE  (RENTRADA.PSD) 


C  THIS  ROUTINE  ESTIMATES  THE  POWER  SPECTRAL  DENSITY  Of 

C  THE  NOISE  IN  AN  IMAGE,  BY  TAKING  THE  MAGNITUDE  (SQUARED) 

C  Of  THE  DISCRETE  SPACE  FOURIER  TRANSfORM  Of  A  (NOISY)  PART 

C  Of  THE  BACKGROUND. 


C  1  CHOOSE  APPROPIATE  REGION,  SET  PARAMETERS  I  AND  J 

C  ACCORDINGLY. 

REAL* 4  RENTRADA(512,5I2),  PSD(32,32),  AR(32,32),  AK32.32), 
BR(32,32),  61(32,32),  ST(32,32),  CT(32,32),  LBR(64) 

C  -  *  *  TAKE  A  PORTION  Of  THE  IMAGE  BACKGROUND  -  *  • 

J-80 

K-80 

DO  L- 1 ,32 
DO  M-1,32 

LL-L+J 

MM=M+K 

AR(L,M>*RENTRADA(LL,MM) 

END  DO 
END  DO 

C  *  »  •  INITIALIZE  IMAGINARY  PART  Of  DATA  *  *  - 

DO  L- 1 ,32 
DO  M-1,32 

AI(L,MH). 

END  DO 
END  DO 


*  *  *  EVALUATE  MAGNITUDE  Of  2-0  FF' T  *  "  * 


CALL  FFTS2(AR,AI,BR,BI,32,32,ST,CT,LBR,64,2,JERR) 

DO  L- 1 ,32 
DO  M-1,32 

PSD(L,MKABS(CMPLX(BR(L,M),BI(L,M))) 


END  DO 


tS3S>EPiai3 


jy*y‘jt,J,v\. 


WV- 


END  DO 


RETURN 


SUBROUTINE  INPUTIM6(INFILENAttE.BYTEIM6.IM6X.lf16Y) 

READ  IN  INPUT  IMAGE  DATA  IN  BYTE  TYPE 
ASSUME  INPUT  DATA  FILE  IS  A  DIRECT  ACCESS  FILE 
BYTE  BYTEIMG(IMGY.IMGX) 

CHARACTER"  40  INFILENAME 


OPEN(UNIT=2  F!LE=INFILENAME,  STATUS='OLD'.  ACCESSOIRECT'  , 
RECORDTYPE=TIXED') 


DO  I-1.IM6Y 

READ(2I)  (BYTElMGtl,  J),  J=1,  IM6X) 

END  DO 

TYPE  » ,  ‘  »  INPUT  IMAGE  HAS  BEEN  READ 

CLOSEC2) 

RETURN 

END 


SUBROUTINE  BYTE_T0JNTE6ER  (BYTEDATA.  INTDATA.  IX.  IY.  IERR) 

THIS  PROGRAM  CHANGE  DATA  (USUALLY  2D  IMAGE  DATA)  IN  THE  BYTE  DATA 
TYPE  INTO  THE  INTEGER  TYPE 


BYTE  BYTEDATA  (IY.  IX) 

INTEGER*  4  INTDATA  (IY.  IX) 

DO  I  =  1.  IY 
DO  J  =  1,  IX 

IF  (BYTEDATA  (I.J)  GE.  -128 
AND.  BYTEDATA  (U)  IT.  0)  THEN 
INTDATAO.J)  =  BYTEDATA  (I.J)  +  256 
ELSE  IF  (BYTEDATA(U)  .GE.O 
.AND.  BYTEDATA(I.J)  IE.  127)  THEN 
INTDATA(U)  =  BYTEDATAO.J) 

ELSE 


I 


TYPE  *  ,'TSDATA  OUT  Of  RANGE!  ROW,  COL,  DATA  =>  ', 

+  I,  J,  BYTEDATA(U) 

I  ERR  =  1 
GO  TO  100 
ENDIF 

END  DO 
END  DO 

TYPE* , '  BYTE  TO  INTEGER  CONVERSION  HAS  BEEN  DONE' 

100  CONTINUE 

RETURN 

END 

SUBROUTINE  INTE6ER_T0_BYTE(INTDATA,  BYTEDATA.  IX.  IY.  IERR) 

C 

C  THIS  PROGRAM  CHANGE  DATA  {USUALLY  2D  IMA6E  DATA)  IN  THE  BYTE  DATA 

C  TYPE  INTO  THE  INTEGER  TYPE 

C 

BYTE  BYTEDATAOX.  IY) 

INTEGER* 4  INTDATAOX,  IY) 

DO  l=1.IX 
DO  J-1.IY 

IF  (INTDATA(U)  IT.  0  .OR.  INTDATA(U)  ,6T.  255)  THEN  TYPE*/* 

DATA  OUTOF  RANGE!  ROW, COL. DATA  => 

+  I,  J,  INTDATAO.J) 

IERR  =  1 
GOTO  100 

ELSE  IF  (INTDATA(U)  ,1E.  127)  THEN 
BYTEDATAO.J)  =  INTDATAO.J) 

ELSE 

BYTEDATA(U)  =  INTDATA(I.J)  -  256 
ENDIF 

END  DO 
END  DO 

TYPE  */  INTEGER  TO  BYTE  CONVERSION  HAS  BEEN  DONE' 

100  CONTINUE 
RETURN 
END 


PROGRAM  LOC-STAT 


C  THIS  PROGRAM  IS  THE  IMPLEMENTATION  Of  THE  IMA6E  ENHANCEMENT  AND 

C  NOISE  FILTERING  ALGORITHM 

C 
C 

BYTE  BYTE  I  MG  (512,512) 

CHARACTER*  1  COEN 

CHARACTER* 40  INFILENAME,  OUTFILENAME 

INTEGER* 2  IMGBUFFC5 1 2,5 1 2),  L AVEfM6(5 1 2,5 1 2).  L0CMEANC5 1 2 ,5 1 2 ). 

-  INT  ARRA  Y(5 1 2,5 1 2) 


C  *  * «  INPUT  FILENAMES  :  INPUT  &  OUTPUT  *  *  * 

TYPE  10 

10  FORMAT  C  INPUT  FILENAME  »>  ',$) 

ACCEPT  15,  INFILENAME 

15  FORMAT  (A40) 

TYPE  20 

20  FORMAT  C  OUTPUT  FILENAME  ■>  ',*) 

ACCEPT  25,  OUTFILENAME 

25  FORMAT  (A40) 

C  ***  INPUT  ORIGINAL  IMAGE  SIZE  *** 

TYPE  50 

50  FORMAT  f  RCW  SIZE  OF  IMAGE  »>  \$) 

ACCEPT  55,  IMGY 

55  FORMAT  (14) 

TYPE  60 

60  FORMAT  f  COLUMN  SIZE  Of  IMAGE  =>  *.$) 

ACCEPT  65,  IMGX 

65  FORMAT  914) 

C  *  *  *  INPUT  WINDOW  SIZE  *  *  “ 

TYPE  66 

66  FORMAT  C  ROW  SIZE  OF  WINDOW  *>  ',$) 

ACCEPT  67,  IWY 

67  FORMAT  (12) 

TYPE  68 

68  FORMAT  ('  COLUMN  SIZE  OF  WINDOW  =>  \$) 

ACCEPT  69,  IWX 

69  FORMAT  (12) 

C  *  ** SELECT  CONTRAST  ENHANCEMENT  OR  NOISE  FILTERING  *  *  * 

TYPE  70 

70  FORMAT  (  CONT RAST/NOISEFIL  ?  (C/N)  =>',$) 


ACCEPT  75,  CE 
75  FORMAT  (A  I ) 


I  ERR  =  0 

<:  88 ‘INPUT  ORIGINAL  IMAGE  888 

CALL  INPUTIM6  (INFILENAME,  BYTEIMG,  IM6X,  IMGY) 

C  8  8  8 CHANGE  IMAGE  FROM  BYTE  TYPE  INTEGER  TYPE 8  8  8 

CALL  BYTE-TOJNTEGER  (  BYTEIM6,  IMGBUFF,  IMGX,  IMGY,  IERR) 

IF  (IERR  EQ.UGOTO  100 

C  888  CALCULATE  LOCAL  MEAN  ARRAY  888 

CALL  LOCAL_MEAN(IMGBUFF,  LOCMEAN,  IMGX,  IMGY.  IWX,  IWY,  IERR) 
IFOERR  .EQ.1  )  GOTO  100 

IF  (CE  ,EQ.  r  .OR,  CE  .EQ.  c )  THEN 

C  8  88  INPUT  K  =  STANDARD  DEVIATION  RATIO  *  8  8 

TYPE  80 

80  FORMAT  C  STAND.  DEV.  RATIO  K=>',  $) 

ACCEPT  85,  K 
85  FORMAT  (14) 

CALL  ENHACONT  (IMGBUFF,  LOCMEAN,  INTARRAY,  IMGX.  IMGY,  K,  IERR) 
ELSE 

C  *  *  *  INPUT  NV  =  NOISE  VARIANCE  *  "  - 

TYPE  90 

90  FORMAT  ('  NOISE  VARIANCE  NV  =>  ',$) 

ACCEPT  95 

95  FORMAT  (14) 

CALL  NOISE-FILT  (  IMGBUFF,  LOCMEAN,  INTARRAY,  NV,  IMGS,  IMGY,  IERR) 
ENDIF 

SCLM1N  «  0.0 
SCLMAX  *  255.0 

C  CALL  SCALE_ADJUST  (INTARRAY,  IMGX,  IMGY,  SCLMAX,  SCLMIN) 

CALL  TRUNCATE  (INTARRAY, IMGX.IMGY) 

C  *  *  *  CHANGE  IMAGE  FROM  INTEGER  TYPE  TO  BYTE  TYPE  *  *  » 

CALL  INTEGER_TO_BYTE  (  INTARRAY,  BYTEIMG,  IMGX,  IMGY.  IERR  ) 

IF  (IERR  .EQ.  1 )  GOTO  100 

C  8  8  8 OUTPUT  THE  PROCESSED  IMAGE  8  8  8 


<  .1"-  sSv'>!  \S  7-  L"  V'-  IT"  i 


CALL  0UTPUTIM6  (  OUTF ILENAME .  BYTEIMG,  IM6X,  inGY  ) 
STOP 

IOO  CONTINUE 

TYPE  •/  !!  FORCED  TO  EXIT  !! ' 

STOP 

END 


SUBROUTINE  ENHACONT  ( IM6BUFF  .LOCMEAN,  INTARRA  Y.IM6X,  IMGY, K.IERR) 


C  THIS  SUBROUTINE  IMPLEMENTS  THE  CONTRAST  ENHANCEMENT  BY  USE 

C  OF  LOCAL  STATISTICS  ALGORITHM. 

C 

C 

INTEGER-2  IMGBUFF(IMGY.IMGX)  .  LOCMEAN(IMGY.IMGX)  , 
INTARRAY(IMGY.IMGX) 

C  - » -  EVALUATE  LOCAL  CONTRAST  AND  MODIFY  IT  WITH  K  -  "  - 

DO  I  -  I  IMGY 
DO  J  ■»  I.IMGX 

INTARRAY(U)  »  K-(IMGBUFF(I,J)  -  LOCMEANO.J)) 

END  DO 
END  DO 

C  -  -  -  ADD  MODIFIED  LOC .  CONTRAST  AND  LOCAL  MEAN  -  -  - 

DO  I  *  1  ,IMGY 
DO  J  -  I  ,IMGX 

INTARRAY  (I.J)  >  INTARRAY  (l,J)  +  LOCMEAN  (l,J) 

END  DO 
END  DO 

C  « -  -  CORRECT  ANY  OVERFLOW  -  »  - 

DO  I  =*  1.IMGY 
DO  J  =  IJMGX 

INTARRAY(U)  =  K-(IMGBUFF(I,J)  -  LOCMEANO.J)  ) 

END  DO 
END  DO 

C  -  -  -  ADD  MODIFIED  LOC .  CONTRAST  AND  LOCAL  MEAN  -  -  - 

DO  I  =  t  .IMGY 
DO  J  =  IJMGX 

INTARRAY  (I.J)  -  INTARRAY  (I.J)  +  LOCMEAN  (I.J) 

END  DO 
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END  DO 


CORRECT  ANY  OVERFLOW  1 


SCLUIN  =  0.0 
SCLMAX  =  255.0 

CALL  SCALE_ADJUST  (INTARRAY,  iM6X.  IMGY,  SCLMAX,  SCLMIN) 


RETURN 

END 


SUBROUTINE  N0ISE_FILT(IM6BUFF,L0CME AN. INTARRAY. NV.IHGX, 
+  IM6Y.IERR) 


THIS  SUBROUTINE  IS  THE  IMPLEMENTATION  OF  THE  NOISE  FILTERING 
ALGORITHM,  CALLED  BY  THE  NOISE  FILTERING  AND  CONTRAST  ENHANCEMENT 
BY  USE  OF  LOCAL  STATISTICS  PROGRAM. 


REFERENCES: 


INTEGER* 2  IMGBUFF(IMGY.IMGX),  LOCMEAN(IMGY.IMGX). 
+  INTARRAY(IMGY.IMGX) 

REAL* 8  AK(  128,1 2d)  .0(128,128) 


*  *  *  EVALUATE  THE  VARIANCE  OF  THE  IMPUT  IMAGE  *  *  * 


TYPE  *  .'CHECK  IMGBUFF  ARRAY’ 
CALL  MCHECK(IMGBUFF) 


TYPE  * , 'CHECK  LOCMEAN  ARRAY  ' 
CALL  MCHECK(LOCMEAN) 

DO  I  =  1.IMGY 
DO  J  =  1  ,IMGX 

IDIF  =  IMGBUFF(U)  -  LOCMEAN(l,J) 

TYPE  16,  IDIF 

FORMAT  ('  DIF  =>  ',110) 

(XI, J)  =  (IDIF*  *2)  -  NV 
TYPE  17,  (XU) 

FORMATC  Q  =>  \F8.2) 

END  DO 
END  DO 


TYPE  *,' CHECK  THE  Q  ARRAY  ' 
CALL  MRCHECK(Q) 

*  *  *  EVALUATE  THF  GAIN  AK  *  *  * 


DO  1=  1.IMGY 
DO  J  =  1,  IMGX 


W'»  -  m  ->  *>  - . 


DEN  =  0(1,  J)  +  NV 
IF  (DEN  .EQ.  0.0)  THEN 
DEN  =  i  .0 

ENDIF 

AK0.J)  =  (0(1, J)  /  DEN) 

END  DO 
END  DO 

C  ***  ESTIMATE  THE  RESTORED  IMA6E  »  *  * 

DO  I  =  1.  IM6Y 
DO  J  =  1 , IMGX 

DIF  =  AK(1,J)** (IMGBUFF(I ,J)  -LOCMEAN0.J)) 
INTARRAY(U)  =  INT(LOCMEAN(l,J)  +  DIF) 

END  DO 
END  DO 

RETURN 

END 


SUBROUTINE  SCALE.ADJUST  (IMGBUFF,  IMGX,  IM6Y,  SCLMAX.  SCLMIN) 

INTEGER- A  IMGBUFF  (IMGX,  IMGY) 

MAXVAL  =  IMGBUFF  (1.1) 

MINVAL  ■  IMGBUFF  (1,1) 

DO  I  =1,  IMGY 
DO  J  =  1 . IMGX 

IF  (IMGBUFFU.J)  .GT.  MAXVAL)  THEN 
MAXVAL  =  IMGBUFF(U) 

ENDIF 

IF  (IMGBUFFd ,J)  .IT.  MINVAL)  THEN 
MINVAL  =  IMGBUFFd, J) 

ENDIF 

END  DO 
END  DO 

C  TYPE  10,  MAXVAL.  MINVAL 

CIO  FORM  ATOM  AX  &MIN  VALUES  =>'  218) 

TYPE  -  ,T1AXVAL=  '.MAXVAL,  T11NVAAL=  '.MINVAL 

INTVAL  =  MAXVAL  -  MINVAL 


SLOPE  =  (SCLMAX  -  SCLMINVREAL  (INTVAL) 


00  1=  I ,  IM6Y 
DO  J  =  1,  IMGX 


ti 


Ik 


I: 


t- 


BS 


EA 


tv 


ft 


IF  (IfiGBUFFd ,J)  .EQ.  MINVAL)  THEN 
IMGBUFF  (l,J)  =  SCLMIN 


ELSE  IF  (IMGBUFF(I.J)  .EQ,  MAXVAL)  THEN 
IMGBUFF  (l,J)  -  SCLMAX 

ELSE 

IDIF  =  IMGBUFF  (l,J)  -  MINVAL 
IMGBUFF(I.J)  =  INT(S)  OPE" IDIF)  +  INT(SCLMIN 

ENDIF 


END  DO 
END  DO 


RETURN 

END 


SUBROUTINE  TRUNCATE0H6BUFF,  IM6X,  IM6Y) 

TRUNCATE  IMAGE  GRAY  VALUES  INTO  0-255 
INTEGER"  A  IMGBUFFOMGY,  IMGX) 


DO  I  =  I,  IMGY 
DO  J  =  1,  IMGX 

IF  (IMGBUFF(U)  .LT.  0)  IMGBUFF(U)  =  0 
IF  (IMGBUFF(I.J)  .GT.  255)  IMGBUFF(1,J)  ■  255 

END  DO 
END  DO 


RETURN 

END 
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PROGRAM  MEDIFILT 


THIS  PROGRAM  IS  THE  IMPLEMENTATION  OF  THE  MEDIAN 
FILTERING  IN  NOISE  REDUCTION 


BYTE  BYTEIMG(5 1 2,  512) 

CHARACTER" 40  INFILENAME,  OUTFILENAME 

INTEGER* 4  IMGBUFF(5 1 2,5 1 2),  IHSTC256),  INTARRAY(512,  512) 


**"  INPUT  FILENAMES  :  INPUT  &  OUTPUT 
TYPE  10 

FORMAT  ('INPUT  FILENAME  =>  \$) 

ACCEPT  15,  INFILENAME 
FORMAT  (A40) 

TYPE  20 

FORMAT  ('OUTPUT  FILENAME  =>  ',$) 

ACCEPT  25,  OUTFILENAME 
FORMAT  (A40) 

" " *  INPUT  ORIGINAL  IMAGE  SIZE  «** 

TYPE  50 

FORMAT  (ROW  SIZE  OF  IMAGE  =>  ',$) 

IMGY=512 
ACCEPT  55,  IMGY 
FORMAT  (14) 

TYPE  60 

FORMAT  C  COLUMN  SIZE  OF  IMAGE  =>  *,  $) 

IMGX  =  512 
ACCEPT  65,  IMGX 
FORMAT  (14) 

TYPE  ",  IMGX,  IMGY 

IERR  =  0 

*  "  *  INPUT  ORIGINAL  IMAGE  "  *  * 

CALL  IMPUTIMG  (INFILENAME, BYTEIMG,  IMGX,  IMGY' 

*  »  "  CHANGE  IMAGE  FROM  BYTE  TYPE  TO  INTEGER  TYPE  *  "  » 

CALL  BYTE_TOJNTEGER  (BYTEIMG,  IMGBUFF,  IMGX,  IMGY,  IERR) 
IF  (IERR  .EQ.  1  )  GO  TO  100 
"  *  *  PERFORM  THE  MEDIAN  FILTERING  "  "  " 

DO  I  =  1 , 20 

CALL  MEDI  (IMGBUFF,  INTARRAY,  IMGX,  IMGY,  5,  5,  IHST,  256) 
DO  IX  =  1,  IMGX 


DO  IY  =  I,  IM6Y 

IMGBUFF(IX.IY)  »  INTARRAYOX,  IY) 

END  DO 
END  DO 

END  DO 

C  **•  CHANGE  IMAGE  FROM  INTEGER  TYPE  TO  BYTE  TYPE  - « * 

C  TYPE  MMGX.  IMGY 

CALL  INTEGER_TO_BYTE  (INTARRAY,  BYTEIMG,  IM6X,  IMGY,  IERR) 
IF  (IERR  .EQ.  1)  GO  TO  IOO 

C  “  *  "  OUTPUT  THE  PROCESSED  IMAGE  *  *  * 


CALL  OUTPUTIMG  (OUTFILENAME,  BYTEIMG,  IMGX,  IMGY) 


STuP 


100  CONTINUE 

TYPE  !  !  FORCED  TO  EXIT  !  !' 
STOP 


■2 


END 
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